CD ROM Paradise Collection 4/CD ROM Paradise Collection 4 1995 Nov.iso/science/cdrift3.zip/TEC2.C

C/C++ Source or Header | 1990-02-22 | 15KB | 291 lines

/* This program is Copyright (c) 1990 David Allen.  It may be freely
   distributed as long as you leave my name and copyright notice on it.
   I'd really like your comments and feedback; send e-mail to
   [e-mail], or send us-mail to David Allen, 10 O'Moore Ave,
   Maynard, MA 01754. */

/* This file contains the function trysplit(), which is called from
   onestep() in tec1.c, and its subfunctions.  One of its subfunctions,
   segment(), is also called from init() in tec1.c.   Trysplit is the
   function in charge of splitting one plate into smaller plates. */


#include "const.h"
#include "var.h"

#define PI 3.14159
#define TWOPI 6.28318
#define TWOMILLIPI 0.00628318

/* RIFTMARK is the temporary indicator placed in the arrays to indicate
   the squares a rift has just appeared.  The function stoprift() puts
   them in, and trysplit() takes them out before anybody can see them. */
#define RIFTMARK -1

/* These are all defined in tec1.c */
extern char m[2][MAXX][MAXY], r[MAXX][MAXY], kid[MAXFRAG];
extern unsigned char t[2][MAXX][MAXY];
extern short karea[MAXFRAG], tarea[MAXPLATE], ids[MAXPLATE], step;
extern struct plate p [MAXPLATE];


trysplit (src) short src; {
   /* Trysplit is called at most once per step in only 40% of the steps.
   It first draws a rift on one of the plates, then it segments the result
   into some number of new plates and some splinters.  If exactly two new
   non-splinter plates are found, new plate structures are allocated, new
   dx and dy values are computed, and the old plate is freed.  If anything
   goes wrong, the rift is erased from the array, returning the array to its
   previous state.  The functions newrift, segment and newplates do most
   of the work. */

   register short i, j, a; short count, old, frag, reg;

   if (old = newrift (src)) if (segment (src, old, &frag, ®)) if (reg > 1) {

      /* Set tarea[i] to areas of the final segmented regions */
      for (i=0; i<=MAXPLATE; i++) tarea[i] = 0;
      for (i=1; i<=frag; i++) tarea[kid[i]] += karea[i];

      /* Give up unless exactly two regions are large enough */
      for (i=1, count=0; i<=reg; i++) if (tarea[i] > MAXSPLINTER) count++;
      if (count == 2) {

         /* Compute new dx,dy, then update m with the ids of the new plates */
         newplates (src, old);
         for (i=0, count=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
            if (a = r[i][j]) m[src][i][j] = ids[kid[a]];
            if (m[src][i][j] == RIFTMARK) {
               m[src][i][j] = 0; t[src][i][j] = 0; } }
         pfree (old); return (0); } }

   /* If execution reaches here, the split operation failed; remove rift */
   for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++)
      if (m[src][i][j] == RIFTMARK) m[src][i][j] = old; }


newrift (src) short src; {
   /* This function randomly picks a center for a new rift, and draws in
   a curving line until the line hits either the coast or another plate.
   If another plate is hit, the rift is invalid and the function returns 0.
   To find a center, the function generates random x,y values until it
   finds one that is at least RIFTDIST squares from any ocean square.  If a
   center is found, a random angle is generated; the rift will pass through
   the center at that angle.  Next, halfrift() is called twice.  Each call
   generates the rift leaving the center in one direction.  If everything
   works out, the function returns the id of the plate the rift is on. */

   short x, y, lx, rx, ty, by, i, j, tries = 0, which; double t;

   /* Generate a random x, y value */
   getctr: if (tries > MAXCTRTRY) return (0);
   x = rnd (XSIZE); y = rnd (YSIZE);

   /* If the location is ocean, try again */
   if (!m[src][x][y]) { tries++; goto getctr; }

   /* Set lx,rx,ty,by to the coordinate values of a box 2*RIFTDIST on a side */
   /* centered on the center.  Clip the values to make sure they are on the */
   /* array.  Loop through the box; if a point is ocean, try another center. */
   lx = (x < RIFTDIST) ? 0 : x - RIFTDIST;
   rx = (x > XSIZE - RIFTDIST - 1) ? XSIZE - 1 : x + RIFTDIST;
   ty = (y < RIFTDIST) ? 0 : y - RIFTDIST;
   by = (y > YSIZE - RIFTDIST - 1) ? YSIZE - 1 : y + RIFTDIST;
   for (i=lx; i<rx; i++) for (j=ty; j<by; j++)
      if (!m[src][i][j]) { tries++; goto getctr; }

   /* Found a good center, on plate `which'.  Put a rift indicator in the */
   /* center.  Generate a random angle t, which is really an integer in the */
   /* range 0-499 multiplied by 2 PI / 1000.  Call halfrift once for each */
   /* half of the rift; t is the initial angle for the first call, and */
   /* t + PI is the initial angle for the second call.  If halfrift() */
   /* returns zero, abort and return 0; otherwise, return the plate id. */
   which = m[src][x][y]; m[src][x][y] = RIFTMARK;
   t = rnd (500) * TWOMILLIPI;
   if (!halfrift (src, x, y, which, t)) return (0);
   t += PI; if (t > TWOPI) t -= TWOPI;
   if (!halfrift (src, x, y, which, t)) return (0);
   return ((int) which); }


halfrift (src, cx, cy, which, t) short src, cx, cy, which; double t; {
   /* Draw a rift from cx,cy on plate `which' at angle t.  At the beginning,
   digitize the angle using Bresenham's algorithm; once in a while thereafter,
   modify the angle randomly and digitize it again.  For each square travelled,
   call stoprift() to see if the rift has left the plate. */

   short ddx, ddy, rdx, rdy, draw, i, a; double dx, dy, adx, ady;

   /* For-loop against SIZE to guard against infinite loops */
   for (i=0; i<XSIZE; i++) {

      /* If first square or 1/6 chance at each step, digitize */
      if (!i || !rnd (BENDEVERY)) {

         /* If not first step, modify angle a little */
         if (i) t = t + (rnd (BENDBY<<1) * TWOMILLIPI) - (BENDBY * TWOMILLIPI);
         if (t > TWOPI) t -= TWOPI; if (t < 0) t += TWOPI;

         /* Compute dx and dy, scaled so that larger is exactly +1.0 or -1.0 */
         dy = sin (t); dx = cos (t); adx = ABS(dx); ady = ABS(dy);
         if (adx > ady) { dy = dy / adx; dx = (dx < 0) ? -1.0: 1.0; }
         else { dx = dx / ady; dy = (dy < 0) ? -1.0: 1.0; }

         /* Convert to integer value and initialize remainder */
         /* for each coordinate to half value */
         ddx = REALSCALE * dx; ddy = REALSCALE * dy;
         rdx = ddx >> 1; rdy = ddy >> 1; }

      /* Main part of loop, draws one square along line.  The basic idea */
      /* of Bresenham's algorithm is that if the slope of the line is less */
      /* than 45 degrees, each time you step one square in X and maybe step */
      /* one square in Y.  If the slope is greater than 45, step one square */
      /* in Y and maybe one square in X.  Here, if the slope is less than 45 */
      /* then ddx == REALSCALE (or -REALSCALE) and the first call to */
      /* stoprift() is guaranteed.  If stoprift returns <0, all is ok; */
      /* if zero, the rift ran into the ocean, so stop now; if positive, the */
      /*  rift ran into another plate, which is a perverse condition and the */
      /* rift must be abandoned.  */
      rdx += ddx; rdy += ddy;
      if (rdx >=  REALSCALE) { cx++; rdx -= REALSCALE; draw = 1; }
      if (rdx <= -REALSCALE) { cx--; rdx += REALSCALE; draw = 1; }
      if (draw == 1) {
         a = stoprift (src, cx, cy, which); if (a >= 0) return (a == 0); }
      if (rdy >=  REALSCALE) { cy++; rdy -= REALSCALE; draw = 2; }
      if (rdy <= -REALSCALE) { cy--; rdy += REALSCALE; draw = 2; }
      if (draw == 2) {
         a = stoprift (src, cx, cy, which); if (a >= 0) return (a == 0); } }
   return (1); }


stoprift (src, x, y, which) short src, x, y, which; {
   /* This function is called once for each square the rift enters.  It
   puts a rift marker into m[src] and decides whether the rift can go on.
   It looks at all four adjacent squares.  If one of them contains ocean
   or another plate, return immediately so that the rift stops (if ocean)
   or aborts (if another plate).  If none of them do, then return ok. */

   register short w, a;

   w = which; p[w].area--; m[src][x][y] = RIFTMARK;
   a = m[src][x][y+1]; if ((a != w) && (a!= RIFTMARK)) return ((int) a);
   a = m[src][x][y-1]; if ((a != w) && (a!= RIFTMARK)) return ((int) a);
   a = m[src][x+1][y]; if ((a != w) && (a!= RIFTMARK)) return ((int) a);
   a = m[src][x-1][y]; if ((a != w) && (a!= RIFTMARK)) return ((int) a);
   return (-1); }


segment (src, match, frag, reg) short src, match, *frag, *reg; {
   /* This routine implements a standard binary-blob segmentation.  It looks
   at the array m[src]; match is the value of the blob, and everything else
   is background.  The result is placed into array r and vectors kid and karea.
   One 8-connected region can be made up of many fragments; each fragment is
   assigned a unique index.  Array r contains the frag indices k, while kid[k]
   is the region frag k belongs to and karea[k] is the area of frag k.
   Variables frag and reg are set on output to the number of fragments and
   regions found during the segmentation.  The private vector kk provides one
   level of indirection for merging fragments; fragment k is merged with
   fragment kk[k] where kk[k] is the smallest frag index in the region. */

   register short i, j, k, k1, k2, k3, l;
   char kk [MAXFRAG];

   /* Initialize all frag areas to zero and every frag to merge with itself */
   for (k=0; k<MAXFRAG; k++) { kk[k] = k; karea[k] = 0; }

   /* Look at every point in the array */
   for (k=0, i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
      /* If too many fragments, give up */
      if (k == MAXFRAG) return (0);

      /* If this square isn't part of the blob, try the next square */
      if (m[src][i][j] != match) { r[i][j] = 0; goto bottom; }

      /* It is part of the blob.  Set k1 to the frag id of the square to */
      /* its left, and set k2 to the frag id of the square above it.  Note */
      /* that because of the for-loop direction, both of these squares have */
      /* already been processed. */
      k1 = i ? kk [r [i-1] [j]] : 0; k2 = j ? kk [r [i] [j-1]] : 0;

      /* If k1 and k2 are both background, start a new fragment */
      if (!k1 && !k2) { r[i][j] = ++k; karea[k]++; goto bottom; }

      /* If k1 and k2 are part of the same frag, add this square to it */
      if (k1 && (k1 == k2)) { r[i][j] = k1; karea[k1]++; goto bottom; }

      /* If k1 and k2 belong to different frags, merge them by finding */
      /* all the frags merged with max(k1,k2) and merging them instead */
      /* with min(k1,k2).  Add k to that fragment as well. */
      if (k1 && k2) {
         if (k2 < k1) { k3 = k1; k1 = k2; k2 = k3; }
         for (l=1; l<=k; l++) if (kk[l] == k2) kk[l] = k1;
         r[i][j] = k1; karea[k1]++; goto bottom; }

      /* Default case is that one of k1,k2 is a fragment and the other is */
      /* background.  Add k to the fragment. */
      k3 = (k1) ? k1 : k2; r[i][j] = k3; karea[k3]++;
   bottom: continue; }

   /* Set up vector kid to map from fragments to regions by using i to count */
   /* unique groups of fragments.  A unique group of fragments is */
   /* characterized by kk[k] == k; otherwise, frag k is merged with some */
   /* other fragment. */
   for (i=0, j=1; j<=k; j++) {
      if (j == kk[j]) kid[j] = ++i;
      else kid[j] = kid [kk [j]]; }

   /* Make sure the id of the background is zero; set up return values */
   kid[0] = 0; *frag = k; *reg = i; return (1); }



newplates (src, old) short src, old; {
   /* Compute new dx and dy values for plates right after fragmentation.  This
   function looks at the rift markers in m[src]; variable old is the index of
   the plate from which the new plates were created.  For each plate adjacent
   to the rift, this function subtracts the number of plate squares to the left
   of the rift from the number to the right; this gives some indication of
   whether the plate should move left or right, and how fast.  The same is done
   for squares above and below the rift.  The results are put into dx[] and
   dy[].  At this point some unscaled movement vector is available for both of
   the new plates.  The vectors are then scaled by the relative sizes of the
   plates.  The idea is that if one plate is much larger than the other, the
   small one should move faster.  New plate structures are allocated for the
   new plates, and the computed dx and dy values are put in them. */

   short dx[MAXPLATE], dy[MAXPLATE];
   register short i, j, a; short totarea=0, maxmag=0; double scale, b;

   for (i=1; i<MAXPLATE; i++) { dx[i] = 0; dy[i] = 0; ids[i] = 0; }

   /* For every point in the array, set a to the region id (kid is the */
   /* lookup table and r contains frag indices); if a is nonzero and */
   /* the rift is adjacent, adjust counters appropriately */
   for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) if (a = kid[r[i][j]]) {
      if ((i-1 > -1)    && (m[src][i-1][j] == RIFTMARK)) (dx[a])++;
      if ((i+1 < XSIZE) && (m[src][i+1][j] == RIFTMARK)) (dx[a])--;
      if ((j-1 > -1)    && (m[src][i][j-1] == RIFTMARK)) (dy[a])++;
      if ((j+1 < XSIZE) && (m[src][i][j+1] == RIFTMARK)) (dy[a])--; }

   /* For those regions larger than splinters (tarea is set up in trysplit), */
   /* allocate a new plate structure and initialize its area; compute the */
   /* magnitude of the dx dy vector and remember the maximum magnitude; also */
   /* record the total area of new regions */
   for (i=1; i<MAXPLATE; i++) if (tarea[i] > MAXSPLINTER) {
      ids[i] = palloc (); p[ids[i]].area = tarea[i];
      totarea += tarea[i];
      a =sqrt ((double) ((dx[i]*dx[i]) + (dy[i]*dy[i])));
      if (a > maxmag) maxmag = a; }

   /* Generate a random speed and predivide so that all speeds computed */
   /* below are less than the random speed. */
   scale = (double) (rnd (SPEEDRNG) + SPEEDBASE) / (maxmag * totarea);

   /* Compute the dx and dy for each new plate; note that the speed the */
   /* plate was moving at before splitting is given by p[old].odx,ody */
   /* but those must be multiplied by MR to get the actual values */
   for (i=1; i<MAXPLATE; i++) if (ids[i]) {
      b = scale * (totarea - tarea[i]);
      p[ids[i]].odx = p[old].odx * MR [p[old].age] + dx[i] * b;
      p[ids[i]].ody = p[old].ody * MR [p[old].age] + dy[i] * b; } }